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The Holstein model of spinless fermions interacting with dispersionless phonons in one dimension 
is studied by a Green's function Monte Carlo technique. The ground state energy, first fermionic 
excited state, density wave correlations, and mean lattice displacement are calculated for lattices of 
up to 16 sites, for one fermion per two sites, i.e., a half-filled band. Results are obtained for values 
of the fermion hopping parameter of t = 0.1a;, w, and 10ui where ui is the phonon frequency. At a 
finite fermion-phonon coupling g there is a transition from a metallic phase to an insulating phase in 
[ which there is charge-density-wave order. Finite size scaling is found to hold in the metallic phase 

OS ■ and is used to extract the coupling dependence of the Luttinger liquid parameters, u p and K p , the 

OS ■ velocity of charge excitations and the correlation exponent, respectively. For free fermions (g — 0) 

and for strong coupling (g 2 3> tu) our results agree well with known analytic results. For t = uj and 
t — Wui our results are inconsistent with the metal-insulator transition being a Kosterlitz-Thouless 
transition. 
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I. INTRODUCTION 



' A wide range of quasi-one-dimensional materials have electronic properties that are dominated by the Peierls or 
OS . charge-density-wave instability caused by the electron-phonon interaction!!! For a half-filled band it is energetically 
favourable for the lattice to dimerize and open an energy gap at the Fermi surface. Although the lattice distortion in- 
creases the lattice energy, opening the electronic energy gap preferentially lowers the total energy for highly anisotropic 
systemsEl. These systems are often modelled by the one-dimensional HolsteinEl or Su-Schrieffer-Heeger (SSH)cl models. 
Most treatments of the Peierls instability treat the phonons in the mean-field or rigid lattice approximation. This 
is questionable in one-dimension and furthermore, in a wide-range of materials the lattice distortion is comparable 
to the zero-point motion of the latticecl. It has recently been- shown that the quantum lattice fluctuations must be 
taken into account to satisfactorily describe opticaLproperties&Q. Several authors have previously considered the role 
of quantum lattice fluctuations for the SSH model! and the Holstein modelErO at half-filling. Voit and Schulz ha*e 
. ,_h , considered the interplay of quantum lattice fluctuations and electron-electron interactions away from half-fillingpJ. 

Recently the Holstein model, also known as the molecular crystal model, has received considerable attention because 
the challenge of high-T c and fullcrcnc superconductors has revealed deficiencies in our understanding of the electron- 
phonon interaction and the competition between superconductivity and charge-cLeHsity-wave instabilities. This has 
motivated, studies of the Holstein model in infinite dimensions!!-], two dimensionsEj, one dimension^ and on just a 



few sites 

We consider the Holstein model in one dimension at half filling and only with spinless fermions, for simplicity. The 
spinless fermions hop along a one-dimensional chain and interact with a phonon mode located on each lattice site. 
The creation operator for a fermion on site i is denoted Cj. The fermions can hop between neighbouring sites with 
amplitude t. In the absence of interactions the phonons all have the same frequency uj, i.e., they are dispersionless. 
The electron-phonon coupling, in units of energy, is g. Phonon position and momentum operators are denoted by 
and pi, respectively. The Hamiltonian for the Holstein model (at half filling) isa 

i i i 

for a system of N lattice sites. This Hamiltonian has particle-hole symmetry since the transformation Cj — > ( — l) z c| , 
qi — > — qi leaves H invariant. This discrete symmetry is broken in the charge-density-wave phase which has the 
electronic order parameter 
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m e=^Y,(- i y <c i*> ( 2 ) 



and the phonon order parameter 



rn. 



!£(-!)<<«> (3) 



which is a measure of the dimerization. 

If phonon creation and annihilation operators are denoted by a\ and Oj, respectively, the Hamiltonian ([l]) can be 
written 

(4°i+i + 4+i c i) - 9^2(c\c t - h + a}) +u)^2,a\ai. (4) 



H 

Thus ground state properties will be determined by two independent parameters, which we will take to be t/w and 
g/uj. It is also useful to define the dimensionless electron-phonon coupling 

A = (5) 

Although for simplicity we confine ourselves to the case of spinless fermions this model is still of physical relevance 
in at least two situations. The first situation concerns strongly correlated electron systems. In the infinite U limit thr; 
Hubbard model is known to map onto the case of spinless-fermionsEfl. This may be realized in the 1^2 TCNQ saltsEj. 
The second situation concerns the spin-Peierls transitioned. Using a Jordan- Wigner transformatiorH this model can 
be mapped onto a XX spin chain in zero field with the Hamiltonian: 

H = -2t £ (SfS? +1 + SfSf +1 ) - g St (a, + a]) + u> £ a\ ai . (6) 

i i i 

It should be pointed out that this is not the standard Hamiltonian used to study the spin-Peierls transition. However, 
it does have the same qualitative features: i.e., a dimerization of the phonons results in a spin singlet ground state 
with an energy gap. _ 

It was recently shownE3 rigorously that for the one-dimensional Holstein model of spinless fermions at half-filling 
there is no long range order for sufficiently small coupling g. Hirsch and Fradkincl studied the Holstein model at 
half-filling using a world-line Monte Carlo technique and a strong coupling expansion. The expansion suggested that 
for spinless fermions quantum lattice fluctuations destroy the dimerized state if the fermion-phonon coupling was 
sufficiently weak and the phonon frequency sufficiently high. The quantum Monte Carlo simulations were performed 
for 0.5lu < t < 3uj and gave a phase diagram qualitatively consistent with the strong coupling expansion. In contrast, 
for fermions with spin their results were consistent with dimerization for finite phonon frequency and all non-zero 
couplings. |— | 

Zheng, Feinberg, and AvignonM used a variational polaron wave function to study the Holstein model at half-filling. 
For spinless fermions the ground state is a charge-density wave for all parameter values. Most of their results were 
consistent with Hirsch and Fradkin. However, they found that for large phonon frequencies (t > 0.3lu) there was 
a first-order phase transition, with a very large jump in the CDW order parameter, between CDW phases when 
g 2 ~ fOw. They point out that this transition may be an artefact of the variational treatment since it is known 
that in small-polaron theory of a single electron a similar two-minimum structure, leading to non-analytic .behaviour 
sometimes referred to as "self-trapping," occurs and is known to be an artefact of the variational teatmentcJ'Eil. 

This paper presents a study of the Holstein model using a Green's Function Monte Carlo technique. Section || 
reviews how the metallic phase shou ld b e a Luttinger liquid and how finite-size scaling can be used to extract the 



Luttinger liquid parameters. Section III briefly summarizes known analytic results of the Holstein model that can be 



used to check and help understand our Monte Carlo results. Section IV contains a detailed description of the Green's 
Function Monte Carlo technique that we use. Our results are presented and interpreted in Section The physical 
picture that emerges from our results is discussed in the final section. 



II. LUTTINGER LIQUIDS AND FINITE-SIZE SCALING 
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A. The Luttinger liquid conjecture 



For weak coupling and high, frequency the system is in a metallic, i.e., gapless phase. According to Haldane's 
"Luttinger liquid" conjectural3'E3 this phase should be in the same universality class as the Tomonaga-Luttinger 
model of interacting spinless fermions. This means the low-energy properties of the metallic phase are completely 
described by an effective Luttinger model with two parameters u p , the velocity of charge excitations or renormalized 
Fermi velocity, and K p , the renormalized effective coupling (stiffness) constant. Important properties of the Luttinger 
model, quite distinct from those of a conventional Fermi liquid, are (i) there are no quasi-particle excitations at the 
Fermi surface and (ii) all correlation functions have non-universal exponents that can be written in terms of the single 
parameter K p . For example, K p determines the singularity of the momentum distribution function close to the Fermi 
surfaceo: 

n(k)~ ~-sign(k-k F )\k-k F \ a (7) 
and of the single-particle density of states 

p(E) ~ \E\ a (8) 

where 

+ £-2). (9) 

For free fermions K p = 1 and a = 0. For attractive (repulsive) interactions K p > 1 (K p < 1). It is remarkable that, 
as explained below, the parameters u p and K p can be determined from numerical calculations on systems of finite 
size. 



B. Finite size scaling 

The ground state energy Eq(N) of a conformally invariant system of N sites is, to leading order in l/iV0, 

(10) 



Eq(N) _ Tru p C 



N " 6iV 2 

where eoo is the ground state energy per site of the infinite system, u p is the velocity of charge excitations, and C 
is the conformal charge. Care must be taken with boundary conditions. We use anti-periodic (periodic) boundary 
conditions for the fermions when there is an even (odd-) number of fermions. This corresponds to periodic boundary 
conditions for the associated spin or bosonic modelsr 1 !. If the system is a Luttinger liquid it belongs to the same 
universality class as the Gaussian model and C = The slope of a plot of Eq(N)/N versus 1/N 2 (compare Figure 
|l|) can then be used to determine u p . 

The energy of the first excited state is, to leading order in 1/N, 

El (N)-E Q (N) = ^ (11) 

where x is the scaling dimension^. A Luttinger liquid has the unusual property that x depends on the coupling 
constants. 

In the presence of particle-hole symmetry x can be related to the correlation exponent K p which determines the 
asymptotic decay of all correlation functions. Let E±i(N) denote the ground state energy of N/2 ± 1 fermions on N 
sites. By particle-hole symmetry E + i(N) = E-i(N). In a general Luttinger liquid of spinless fermionsta with charge 
density n the compressibility k is given by 



1 _ d 2 e oc (n) _ mi p 
i 2 k dn 2 K p 



(12) 



Since particle-hole symmetry implies de ^°^ = it follows that 



£ ±1 (A)^ (A) + i(l) 2 A%^ (13) 
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which with ( |l2| ) implies that to leading order in 1/N 



E ±1 (N)-E (N) = ^^ (14) 

This is identical to ( pi] ) with K p — 1/Ax. Hence, if u p is known a plot of the energy gap versus 1/N (compare Figure 
||) can be used to determine K p . 

III. ANALYTIC RESULTS 

Certain limits of the Holstein model for which analytic results can be obtained are now briefly reviewed. These 
results will be compared to the appropriate numerical results. 

A. Localized fermions (t = 0) 

The fermions cannot move between sites and the Hamiltonian reduces to N independent Hamiltonians. The 
Hamiltonian for the ith site is 

Hi = + \Mu 2 qf - (m - \)g{2Muj) 1 /\ - (15) 

where rii = cjcj is the fermion occupation at site i. The presence or absence of a fermion shifts the equilibrium 
position of the oscillator to +q e or — g e , respectively, where 

{ M \ * 

qe=g \2^) = < *( 2rl «- 1 )> ( 16 ) 

This Hamiltonian can be diagonalized by the Lang-Firsov transformation^: c, — > c, exp(q e (a\ — Oj)), — > — (2n^ — 
l)g e . The mean square lattice displacement is 



The ground state energy per site is 



4:Uj' 



(18) 



B. Small Polarons (g 2 > tuj) 

This corresponds to the case of a narrow band of small polaronsi. The intcrsite hopping represents a small pertuXi. 
bation on the situation considered in the previous section. Hirsch and Fradkin[!l derived an effective HamiltonianLj 
involving the new fermion operator C[ that creates a fermion at site i and changes the oscillator ground state from 
one centered at —q e to one centered at +q e . 

Heff - -iv£ (c?0+i + C} +I a) (cJCi - (cj +1 C i+1 ^ (19) 

The first term is the polaron binding energy (compare equation (|l8|)) and dominates the ground state energy (Figure 
^) . The second term describes hopping between neighbouring sites with the bandwidth reduced by the overlap of the 
oscillator ground state centered at —q e and +q e : 

J = texp(-(£) ). (20) 



4 



The third term describes the second order process (of order t 2 /uj) where a fermion hops to a neighbouring site and 
back again. This term is repulsive because this process is not possible if the neighbouring site is occupied. 

V=— dx- (21) 

w Jo x 

There is an additional term, of second order in t 2 /w, involving next nearest neighbour inteactions but it is smaller by 
a factor of order 1/A and so the effective Hamiltonian should accurately describe the physics in the strong coupling 
limit, A|-f^ 1 or g 2 ^> tuj. Note that the limit u> — > 00 corresponds to free fermions, as pointed out by Hirsch and 
FradkirM 

The effective Hamiltonian (|l9| ) is, after a Jordan-Wigner transformation, of the same form as that of the exactly 
soluble antiferromagnetic XXZ quantum spin chained. For convenience we now briefly summarize some of the known 
results for this model. It can be exactly solved by Bethe ansatzErEj. The system is metallic for V < 2 J, i.e., there is 
no energy gap or long range order and so it is a Luttinger liquid. It is in an insulating charge-density-wave state for 
V > 2 J. The metal-iasulator transition is an infinite order, i.e. Kosterlitz-Thouless, transition and has been discussed 
in detail by ShankaiO. 

Define a new variable fi by 

V 

cos^=— (22) 

where < /1 < ir/2. As V increases from to the transition at V — 2J, /1 decreases from ir/2 to zero. The velocity 
of charge excitations is given byc3 

u p = ttJ . (23) 

A* 

As V increases from to the transition at V = 2J, the velocity increases from 2 J to ir J. However, as the coupling g 
increases J rapidly decreases and so u p rapidly decreases (compare Figure (a)). The Luttinger liquid exponent K p 



As V increases from to the transition at V — 2 J, K p decreases from 1 to 1/2 (compare Figu*a-0 (b)). The value 
K p = 1/2 is a universal feature of a Kosterlitz-Thouless transition for one-dimensional fermionsT 
For V larger than 2 J define a new variable 7 by 



5 3ll 



The charge-density-wave order parameter (H) is 



V 

cosh 7=—. (25) 



in- =1 tanh 2 (m 7 ). (26) 



The coupling dependence is shown in Figure^ (a) for t = O.Ioj. The metal-insulator transition is Kosterlitz-Thouless 
although it does not appear so on the scale shown. The energy gap in the insulating phase is 

A = 2Jsinh 7 V \ / - (27) 
cosh(m7) 

m=— 00 

and turns out to be very small (compare Figure |5| (b)). 

C. Free fermions (g = 0) 

The fermion states are plane waves with energy dispersion 

E(k) = -2*cos(fc). (28) 

These states are occupied for |fc| < kp = ir/2. Near the Fermi surface at k = ±kp we have E(k) — ±2t(k =p ^f) and 
so the Fermi velocity is Vp = 2t. The ground state energy per site is 

„ r' 2 dk 2t 
e 00 = -2t/ — cos(fc) = . (29) 

J— it/2 27r T 
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D. Adiabatic or mean- field limit (lo < ie~ 1/A ) 



It is assumed that the fluctuations of the lattice about its dimerized value can be neglected and the quantum 
operator qi in the Hamiltonian (|l|) is replaced by its mean value: qi — >< qi >= (—l) l m p . The fermionic Hamiltonian 
can then be diagonalized by a Bogoliubov transformation and the fermionic energies are 

E{k) =±((2tcos(fc)) 2 + A 2 )* (30) 

where A = g(2AIcj) 1 ^ 2 m p is the energy gap at the Fermi surface due to the dimerization. A is then treated as a 
variational parameter and the total energy of the system is minimized to give the self consistent equation 

r /2 i 

1 = At / dk r . (31) 

J-n/2 ((2i cos(fc)) 2 + A 2 ) 5 

The system is dimerized for all coupling strengths and for weak coupling (A < 1) the energy gap is 

A = 8iexp( — -]. (32) 



v A 

The charge-density-wave order parameter is 



— expl-i). (33) 



2lT\t 7rA V ^ 

The corrections to the mean-field equation (|3l|), to next order in uj\/A, were recently calculatedll 



IV. THE GREEN'S FUNCTION MONTE CARLO METHOD 



At first we tried simulating the model using a discrete basis of free phonon eigenstates on each site and employing 
a "stochastic truncation"l2a technique appropriate to this basis. This method gave accurate results for small coupling 
g, but not at or beyond the metal-insulator transition. In this region the staggered displacement m p becomes large, 
corresponding to the presence of highly excited states in the free phonon eigenstate basis. It was thus found more 
appropriate to use a continuous "position space" basis with variables {qi}, and use a different Monte Carlo technique 
as described below. 



A. Ground state energy 



To simulate-the model, we use a Green's Function Monte Carlo (GFMC) method, as develoned by Kalos and 
collaboratorsLZlH, and applied to lattice gauge theory by Chin, Negele and KooninE£l and otherscljllj. Let us review 
the method briefly. 

The Hamiltonian for the Holstein model ([j]) can be rescaled to the dimensionless form 

H = ( c k+i + 4+iCi) + Y,Pi + 9? - 9^(n t - \) (34) 

i i 

where t = 2t/co, g = 2y2g/uj. The imaginary-time Schrodinger equation for the system reads 

- -^\$(t) >= (H - E t )\$(t) > (35) 

where Et is a trial energy, representing a constant shift in the zero of energy. Evolving this equation for time Ar 
yields 

|$(r + At) >= exp (At(E t - H)) |$(r) > . (36) 
At large times r the ground state will dominate: 



G 



|$(t) >~ c exp(-(_B - E t )t) |$ > as r — > oo (37) 

where |$o > is the (time-independent) ground state of -ff with energy E . 
We shall work in a position-space representation, where the wave function 

$({< Z ,n},r)=<{< Z ,n}|<f(r) > (38) 

and \{q, n} > represents an eigenstate of the positions {qi} and fermion occupation numbers {rii} at each site. In this 
representation, 

H = H Q + H 1 (39) 

where 

Hi = -tJ2 + 4f iCi) (4°) 

i 

is the fermion hopping term, and 

^-eS+^-d (41) 

i % 

with 

% i 

as the "potential" term. 

The evolution equation (^) now has the form of a diffusion equation in configuration space. It is assumed that the 
ground-state wave function can be chosen positive everywhere, and it is simulated by the density distribution of an 



ensemble of random walkers in configuration space, whose time evolution mimics that of equation (36). 

To obtain good accuracy, one needs to introduce some variational guidance, which can be done as follows. Let 
> be a trial vector, e.g., some variational approximation to the true ground-state eigenvector with wave function: 

V T ({q,n}) =< {q,n}\V T > . (43) 

Then carry out a similarity transformation 

|$(r) >-> |$'(r) >= W t |$(t) > (44) 

H -> H' = ^tH^ 1 (45) 

where the transformation matrix is diagonal in {q, n} space, with diagonal entries ^^({q, n}). The modified 
evolution equation will be 

|$'(r + At) >= exp (At(E t ~ H')) |$'(r) > (46) 

Let us now separate the fermion hopping term from the rest of the Hamiltonian, and write for small At 

exp {At{E t ~ H')) ~ exp (At{E t - H' Q )) [1 - AtH[] + 0(Ar 2 ) (47) 

(All our calculations from here on will only be accurate to O(Ar)). 
Now Hq transforms to 

i 1 

= " + ^ (-5S-) ^ + {-afr ) ] + v{{q > n}) 

= V^HoVt + + 2i P* (^t 1 ^ 1 -)} (48) 
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as shown by Chin, Negele, and KooninEa where the operator pi = — i-J^ acts on everything to the right of it as 
usual. Then the matrix clement between position eigenstates corresponding to the time-step At at iteration m can 
be writteno 



< {q, | exp (Ar(E T - H' )) \{q, n}^ > 



i 7 j I Mr \ \ J- U 



(m+1) - ^ (m) - 2Ar*: 



0(At 2 ) 



(49) 



Representing the wave function by a distribution of random walkers in position space, the Monte Carlo simulation 
proceeds as follows. Each iteration corresponds to a time step At, and involves a sweep through each site in turn. 
The first term in the exponential (E9J) is simulated by a displacement of each position variable 



Aq, = 2At^t — ± +X 



(50) 



where x is randomly chosen from a Gaussian distribution with standard deviation V2Ar. The first term in ( pp| ) is 
the "drift" term, and the second is the "diffusion" term. The second term in the exponential (^) is simulated by 
multiplying the "weight" of each walker by an equivalent amount. 
We also need to simulate the fermion hopping term: 



< {q, I [1 - AtH[] \{q, n}^ > 

* T ({g,n}(" l+1 )) 



* T ({g,«} (m) ) 



<{g,n}( m+1 )|[l+fAT^(clc 4+1 



\\{q, n ) 



(m) 



> 



(51) 



The factor in front produces a "reweighting" of the walkers in the ensemble; while the hopping term itself produces 
new configurations on walkers with different fermion occupation numbers. 

At the end of each iteration, the trial energy Et is adjusted to compensate for any change in the total weight 
of all walkers in the ensemble; and a "branching" process is carried out, so that walkers with weight greater than 
(say) 2 are split into two new walkers, while any two walkers with weight less than (say) 0.5 .are combined into one; 
chosen randomly according to weight from the originals. This procedure of "Runge smoothing" c3 maximises statistical 
accuracy by keeping the weights of all walkers within fixed bounds, while minimizing any fluctuations in the total 
weight due to the branching process. When equilibrium is reached after many sweeps through the lattice, the average 
value of the trial energy Et will give an estimate of the ground-state energy Eq, from equation ([37]); and the density of 
the ensemble in configuration space will be proportional to $o , J / t- Various corrections due to the finite time interval 
At have been ignored in this discussion, and the limit At — > must be taken in some fashion to eliminate such 
corrections. 

As a trial wave function, we choose a Gaussian, displaced by an amount go at each site depending whether the site 
is occupied or unoccupied: 



^t(U, n}) = exp 



■ ^2(q t - 2q (m - ]-)f 



where c and qo are variational parameters. Then the local "trial energy" 

E L ({q, n}) = V^HoVt = - 9Qi(ni - h - $>c 2 fe - 2q (n i - \)) 2 - 2c] 



(52) 



(53) 



and the "drift" term is 



2* 



dqi 



-4c(q t - 2q Q (n l - -)) 



while the "reweighting factor" in equation ( |5l|) is 

* T ({<?, n}(™ +1) ) 
* T ({<?,fi} (m) ) 



exp 



(m+l) „( m )^ 



(54) 



(55) 
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If the choice of trial function is a good one, and Et is adjusted to approximately equal E$ , then we will have 

E L ~ E T ^ E (56) 

so that the weight of each walker changes very little at each time step, according to equation (|l|), so that the 
fluctuations in the weights are small, and consequently the accuracy of the calculation is maximised. 



B. Expectation Values 



pectation values can also be measured, using a "secondary amplitude" technique discussed by 
I. Let < Q >o be the expectation value to be measured, where we assume the operator Q is 



(57) 



Ground-stat* 
Hamer et alS3 

diagonal in the {q, n} representation. Use Q as a perturbation to the Hamiltonian: 

H' = H + xQ. 



Let Eq(x) denote the ground state expectation value of this Hamiltonian. By the Hellmann-Feynman theorem, the 
required expectation value is given by 



< o > - dE '° 

< H >0— — 

ax 



a:=0 



Taylor expand the eigenvector and eigenvalue 

|$(r,x) >= |$°(r) > +x\$ 1 (t) > +0{x 2 ) 



(58) 



(59) 



E' (x) = E + xE 1 + 0(x 2 ) (60) 

substitute in the evolution equation ( |36|) (ignoring any variational guidance for the time being), and equate powers 
of x to obtain: 

|$°(r + At) >= cxp (At(E t - H)) |$°(r) > (61) 

and 

|$ x (r + At) >= exp (At(E t ~ H)) ^(t) > +At(E 1 - Q)|$°(r) > . (62) 

Equation ( |6l|) is just the original evolution equation (|3^) for the unperturbed system. Equation ( |62| ) is an evolution 
equation of similar structure for the first-order wave function l^ 1 >. It is simulated by giving a "secondary" weight to 
each walker in the ensemble, and evolving it according to (^); while a secondary trial energy E' T is used to estimate 
E 1 , and is adjusted after each iteration to compensate for any change in the total of all secondary weights. At 
equilibrium, the average value of E' T gives an estimate of E 1 , which is equivalent to < Q >o by equation (p8[). 
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V. RESULTS AND DISCUSSION 

GFMC runs were performed for a range of differ- 
ent couplings g/tu at hopping parameter values t — 
0.1cl>,cl> and lOoi, for lattice sizes of 2, 4, 6, 8, and 16 
sites. In each case, the variational parameters c and qo 
(compare equation (^2|)) were adjusted to their optimum 
values by a series of trial runs. Production runs typi- 
cally employed an ensemble of 2000 walkers for 20,000 
iterations. The first 2000 iterations were discarded to 
allow for equilibrium, and the remainder were averaged 
over blocks of 1024 iterations before estimating the error 
to minimize correlation effects. Calculations were per- 
formed on a cluster of six HP 735 workstations. A typical 
run for 16 sites took 1-2 hours of CPU time. Two differ- 
ent time steps were used in each case, namely At = 0.005 
and 0.01 at t = O.lu, At = 0.0005 and 0.001 at t = u, 
and At = 0.001 and 0.002 at t = lOw. The results were 
then linearly extrapolated to At = 0. 

The quantities measured were the ground-state en- 
ergy (in the half-filled sector), Eq(N), the energy gap 
(to the "one hole" sector with one fewer fermions), 
E^i(N) — Eq(N), and ground-state expectation values 
for the mean displacement < ft (2^ — 1) >, the mean 
square displacement < qf >, and two correlated fermion 
expectation values, < n-i«.i+jv/2 > an d < n i n i _ 1+N ^ 2 >, 
where N is the lattice size. The difference between these 
last two values provides an estimate of the amount of 
"staggering," or dimerization, in the fermion occupation 
numbers. (Compare equation ( |6^ ) below). A sample of 
results is shown in Table |[ 

The charge velocity u p was extracted from a finite size 
scaling plot of the ground state energy per site Eq(N)/N 
versus 1/A 2 (compare Figure |l]). According to equation 
( |l0| ) this should be a straight line for large N. To allow 
for the small curvature of our plots, because we used 
only moderately large system sizes (N = 2, 4, 6, 8, and 
16 sites), the data was fitted to 



Ep{N) 
N 



6N 2 



a 



(63) 



The correction-to-scaling term 0(N Z^X-m&tches that 
predicted to hold for the XXZ modelEIrca, at least for 
weak interactions (i.e., small /x). For stronger inter- 
actions the exponent is interaction dependent. At the 
metal-insulator transition the correction-to-scaling term 
is 0((N\nN)~ 2 ). However, we found that using such a 
form did not improve the quality of the least square fits 
near the transition. 

For free fermions (g — 0) the values of and u p 
extracted from the fits were found to agree well with the 
known analytic results — — Itjix and u p — 2t (Table 
|J). For t = O.lw the dependence of the ground state 
energy on the coupling is in good agreement with small 
polaron theory (Figure |J). 

The energy gap A of the infinite system was extracted 
from finite size scaling plots of the hole energy E_i(N) — 



Eq(N) versus 1/N (compare Figure g). To allow for the 
corrections to scaling this was fitted to 



E^i(N) — Eq(N) = A + 



1 



2K P N 



w (64) 



Again, the higher order term 0(N~ 3 ) was chosen to be 
consistent with known results for free fermions and the 
XXZ model with weak interactions. To extract K p we 
need to use the value of u p extracted earlier. (Strictly 
speaking equation ( |64j ) is only valid when A = but wc 
use A as a parameter in our fits to check that we are in 
the critical regime. Also, the derivation of equation d&jj ) 
requires particle-hole symmetry, i.e., E_i(N) = E + i(N). 
We checked for several parameter values that the Monte 
Carlo results were consistent with this.) 

Figure ^ (a) shows the dependence of the charge ve- 
locity u p on the fermion-phonon coupling g for t = 0.1a;. 
The results are in good agreement with equations fl2C| ) 
and (|23| ) (solid line in Figure §(a)). The charge veloc- 
ity is significantly reduced by polaronic band narrowing. 
The correlation exponent K p is shown in Figure [| (b) 
as a function of the fermion phonon coupling. The de- 
pendence of K p on the coupling is consistent with the 
metallic phase being a Luttinger liquid. The fact K p < 1 
indicates repulsive interactions in the Luttinger liquid. 
Kp is not plotted for g > 1.5w because the relative error 
is very large. This is because its determination depends 
on the value of u p which has a very large relative error for 
g > 1.5a; (see Figure [| (a)). For t = O.lu;, equations ( p0| ) 
and (|l]) , for the small polaronic model together with the 
criterion V — 2 J can be used to determine that the tran- 
sition from the Luttinger liquid to the insulating phase 
occurs when g = 2.075a;. 

The charge-density- wave order parameter m e , defined 
by (3), must be zero for any finite size system. However, 
in the dimerized phase we also have for j large 



< mni +j >= - + (~l) J m e 



(65) 



and so 



m e = 2 1 < n i n i+N/2 > 



<n i n i „ 1+N/2 > \. (66) 



This equation was used to determine m\ from the results 
for iV = 16 sites. 

Figure || (a) shows the coupling constant dependence 
of m\ for t — O.lu;. The quantum Monte Carlo data 
suggests there is a transition near g = 1.8a;. This is 
consistent with the small polaron theory prediction of 
g = 2.075a; since the latter theory is only valid to order 
1/A ~ Tituj/g 2 , i.e., about 10%. Figure || (b) shows the 
energy gap as a function of coupling. It is not possible to 
detect the transition in the energy gap data. Small po- 
laron theory predicts an energy gap smaller than typical 
uncertainties in the Monte Carlo data. 

Figure H shows the coupling strength dependence of the 
mean lattice displacement < qi(2n,i — 1) > and the mean 



10 



square lattice displacement < qf > for t = O.lui and a 
system of 16 sites. The results are very close to those an- 
ticipated for localized fermions (compare equations Jl6| ) 
and (fl7l)). For t = uj and t — IOuj the mean lattice dis- 
placement was also non-zero, i.e., the ground state was 
polaronic for all couplings, although the magnitude of 
the displacement decreased significantly with increasing 
t. Similar trends are seen for the Holstein model on two 
sitesca. 

For t = uj the charge velocity is again reduced by po- 
laronic effects (Figure fj] (a)) but not by as much as for 
t = Q.Iuj. The interactions in the Luttinger liquid are 
now attractive (K p > 1). Both the order parameter and 
the energy gap show a transition to the insulating phase 
near g — 1.7 'uj (Figure ||). Clearly our results are incon- 
sistent with the universal value .K„ = 1/2 expected for a 
Kosterlitz-Thouless transitionoa. 

For t = 10a; the charge velocity decreases by less than 
ten per cent with increasing g/uj (Figure ^J). This is in 
contrast to the cases of t = O.Iuj and t = uj for which 
the charge velocity decreases by about an order of mag- 
nitude. The correlation function exponent K p increases 
by about fifty per cent. As for t — uj, K p ^ 1/2 at the 
transition and so the transition cannot be a Kosterlitz- 
Thouless transition. The order parameter m e and energy 
gap A become nonzero about g = 3.5a; (Figure |To|). This 
is quite .different from what is anticipated by Hirsch and 
Fradkircl. They performed simulations from t — 0.5a; up 
to t = 3.1a;. They found a smooth decrease in the critical 
value of A c = g 2 /irujt with increasing t/uj, and anticipated 
a smooth crossover to A c = for uj = 0. Extrapolating 
their results to t — \Quj gives A c ~ 0.01 and g c ~ 0.6a; 
compared to our value of g c ~ 3.5a;. Note that the ratio 
of the energy gap to its mean field value (Figure [l0| (b)) 
is much smaller than the ratio of the charge-density-wave 
order parameter to its mean field value. This is consis- 
tent with work showing that the zero point motion of the 
lattice can reduce the magnitude of the order parameter 
by a small amount but produce a-,substantial subgap tail 
in the fermionic density of statesu. (For example, results 
on the continuum version of-lhe SSH model shown in 
Figures 1 and 3 of ReferenceO show that for one set of 
parameter values the energy gap can be about 60 % of 
the mean-field value while the order parameter is only 
reduced by about 5 %). 

The phase boundary as a function of t/uj and g/to be- 
tween the metallic and insulating phases is shown in Fig- 
The solid curve is the boundary predicted by 
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IIIB 



ure 

small polaron theory and the XXZ model (Section 
This curve is only shown for t < uj since this model is 
only valid in the strong coupling limit (t <C g 2 /uj). The 
crosses are the boundary points deduced from Figures 
|8|, and [h]. It should be stressed that there is some ambi- 
guity in determining the phase boundary. According to 
mean-field theory the transition occurs at g = but the 
solid curves in Figures || and [l(] suggest that the transi- 
tion is actually only detectable at g ~ 0.6a; and g ~ 2a;, 
respectively. On the other hand, for t = O.luj small po- 



laron theory and the XXZ model predict a Kosterlitz- 
Thouless transition at g = 2.075a; and the solid curve 
in Figure || (a) shows there is very little ambiguity as- 
sociated with this transition point. For comparison the 
boundary points found by Hirsch and Fradkincl (Figure 
11 in their paper) are also shown. For t = uj there is a 
discrepancy between our results and theirs: they observe 
the transition at smaller coupling than we do. We have 
no explanation for this discrepancy. 



VI. CONCLUSIONS 

We have shown that the Green's function Monte Carlo 
technique can be successfully used to investigate a one 
dimensional fermion-phonon model. As far as we are 
aware this is the first application of this technique to this 
important class of models. The results were of sufficiently 
high precision that a finite size scaling analysis of the 
results could be performed. For the case of free fermions 
(g = 0) and the strong coupling limit (g 2 >> tuj) our 
results agree with known analytic results. 

Our results are consistent with the following physical 
picture of the Holstein model of spinless fermions at half- 
filling. For sufficiently weak coupling the system is in 
a metallic, i.e., gapless phase, with the properties of a 
Luttinger liquid, i.e., the exponents associated with the 
decay of correlation functions depend on the coupling 
strength. The fermions are polaronic, i.e., there is a fi- 
nite phonon displacement q e associated with the occupa- 
tion of a site by a fermion and the velocity of excitations, 
Up, is reduced below the free electron value 2t. As the 
coupling g increases and t/uj decreases q e increases and 
Up (which is a measure of the polaronic band width) de- 
creases. Qualitatively spiilar behaviour is seen for the 
two-site Holstein modeled. In the anti-adiabatic limit 
(t<w) the effective interaction between polarons is re- 
pulsive and for strong coupling the Holstein model can 
be mapped onto the XXZ antiferromagnetic spin chain 
(Section |III B ). However, as t/uj increases the effective 
interaction becomes attractive. This is indicated by a 
change in the value of the stiffness constant K p from val- 
ues less than one to values larger than one. When the 
fermion-phonon coupling is sufficiently large the system 
undergoes a transition to an insulating phase, i.e., an 
energy gap opens at the Fermi surface. There is long- 
range charge-density-wave order and a dimerization of 
the phonons in this phase. Our results for t = uj and 
t = \Quj are inconsistent with the metal-insulator transi- 
tion being infinite order. On the other hand we do not 
see any evidence of the first-order, transition suggested 
by Zhenft. Feinberg, and AvignonEil and by Wu, Huang, 
and SunE3 for certain parameter values. 

This work suggests several possible future investiga- 
tions which we plan to pursue: (a) The adiabatic re- 
gion of the phase diagram (t uj), in which we found 
a larger region of the metallic phase than anticipated by 
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Hirsch and Fradkin, needs to be investigated further, (b) 
The relative importance of superconducting and charge- 
density-wave correlations should be investigated in the 
region of the metallic phase for which the effective inter- 
actions are attractive, (c) Alternative variational wave- 
functions, such as the double Gaussian proposed by Shore 
and Sandem, could be used instead of the single Gaus- 
sian (|52| ) used for the variational guidance. Finally, we 
plan to use this method to investigate the Holstein model 
with spin, and away from half-filling, as well as the Su- 
Schricffcr-Heeger model, and the spin-Peierls problem. 
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FIG. 1. Finite-size scaling of the ground state energy 
Eo(N) for different values of the fermion-phonon coupling 
g. The data shown are for TV = 4, 6, 8, and 16 lattice 
sites. All data are for a fermion hopping parameter t equal 
to the phonon frequency u and for a half-filled band (i.e., one 
fermion per two sites). All energies are in units of u. If the 
system is critical for a particular g value then the data for N 
large should lie on a straight line (see equation ([H^)). The 
lines are least square fits to a parabola (see text). The errors 
in the Monte Carlo data are smaller than the symbol sizes. 




FIG. 2. Finite-size scaling of the hole energy for different 
values of the fermion-phonon coupling g. The data shown are 
for TV = 4, 6, 8, and 16 lattice sites. Eq(N) is the ground 
state energy of a system of N/2 fermions and E-i(N) is the 
ground state energy of a system of N/2 — 1 fermions. All data 
are for a fermion hopping parameter t equal to the phonon 
frequency u). All energies are in units of ui. If the system is 
critical then for that g value the data for large TV should lie 
on a straight line through the origin (see equation (^) ) . The 
lines are least square fits to a cubic (see text). 



FIG. 3. Dependence of the ground state energy per site 
too on the fermion-phonon coupling g for t = O.Icj. The solid 
lines are the predictions of the small polaron model (Section 
[II B). The error bars are smaller than the symbol size. All 
energies are in units of ui. (a) eaa is deduced from the inter- 
cept of the finite-size scaling plot of the ground state energy 
(compare Figure [l]). The dashed line is the polaron binding 
energy —g 2 /Auj and clearly represents almost all of the ground 
state energy, (b) The polaron binding energy has been sub- 
tracted from the ground state energy to make a more accurate 
comparison with small polaron theory. 




FIG. 4. Dependence of Luttinger liquid parameters on 
the fermion-phonon coupling g for t = O.lu. The solid lines 



are the predictions of the small polaron model (Section III B ) . 
(a) The velocity of charge excitations u p is deduced from the 
slope of the finite-size scaling plot of the ground state energy 
(compare Figure Q) and is normalized by the free fermion 
value 2t. The decrease of u p with increasing g is due to the 
narrowing of the bandwidth by polaronic effects, (b) The 
correlation function exponent K p is deduced from the ratio of 
the slopes of the finite-size scaling plots in Figures |l] and ^ (see 
equation (|TT|) ) . The error bars are based on the uncertainties 
in the least-squares fits to the finite size scaling data (compare 
figures [l] and ^) . 
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FIG. 5. Dependence on the fermion-phonon coupling g of 
(a) the square of the charge-density-wave order parameter m e 
and (b) the energy gap A. The solid lines are the predictions 
of the small polaron model (Section [II B). It predicts an in- 



finite order transition at g — 2.075aj. mi was deduced from 
equation (^6|) for a system of 16 sites. The energy gap was de- 
duced from the N — oo extrapolation of the finite size scaling 
plot of the hole energy (compare Figure ^). Both the order 
parameter and the energy gap become non-zero for g > 1.8 
marking the transition into the insulating phase. 
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FIG. 6. Dependence on the fermion-phonon coupling g of 
the mean lattice displacement < qi(2n; — 1) > and the mean 
square lattice displacement < q\ > for t — O.lw and a system 
of 16 sites. The displacements are in units of (Ma>) _1/2 . The 
solid lines are the predictions for localized fermions (Section 




g/<3 

FIG. 7. Same as Figure | but with t = u. The fact 
that K p depends on g and is larger than one is consistent 
with the metallic phase being a Luttinger liquid with attrac- 
tive interactions. If the metal-insulator transition was Koster- 
litz-Thouless K p would equal 0.5 at the transition. 
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FIG. 8. Same as Figure [B| but with t = oj. The solid curves 
are the predictions of mean field theory (Section HID). 
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FIG. 9. Same as Figure ^ but with t = lOw. Note that the 
vertical scale is expanded compared to Figures ^ and ^ Since 
K p 7^ 0.5 near the transition the transition is not Koster- 
litz-Thouless. 
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FIG. 10. Same as Figure || but with t = 10w. Th e soli d 
curves are the predictions of mean field theory (Section [IIP ). 
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FIG. 11. Phase diagram showing the boundary between 
the metallic and insulating phase. The solid curve is the pre- 
diction of small polaron theory and the XXZ model (Section 
MB). The crosses are the results of thisistudy and the solid 
squares the results of Hirsch and FradkirH 
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TABLE I. Monte Carlo results for different quantities for t = u> and g = 1.5a; and for various system sizes. The energies are 



in units of u> and displacements in units of (Mlj) 



N 


E B (N)/N 


£_i(JV) - E„(N) 


< qi(2m - 1) > 


< 1i > 


< ™;™i+M/2 > 


< n i n i _ 1 + M/2 > 


2 


-1.143 ± 0.001 


1.158 ± 0.004 


0.313 ± 0.005 


0.748 ± 0.008 


0.000 ± 0.000 


0.500 ± 0.000 


4 


-0.895 ± 0.001 


0.416 ± 0.004 


0.446 ± 0.001 


0.864 ± 0.004 


0.260 ± 0.001 


0.120 ± 0.0004 


6 


-0.868 ± 0.001 


0.268 ± 0.009 


0.470 ± 0.005 


0.883 ± 0.007 


0.217 ± 0.001 


0.244 ± 0.001 


8 


-0.861 ± 0.002 


0.200 ± 0.015 


0.484 ± 0.002 


0.902 ± 0.003 


0.240 ± 0.002 


0.229 ± 0.001 


16 


-0.854 ± 0.001 


0.064 ± 0.027 


0.488 ± 0.004 


0.904 ± 0.005 


0.247 ± 0.002 


0.246 ± 0.002 



TABLE II. Comparison of Monte Carlo results with known results for free fermions. The ground state energy per site of 
the infinite system, e^, and the velocity of charge excitations, u p , are normalised to their free fermion values. The correlation 



exponent K p is one for free fermions. 



t 


2t 


Up 

It 


Kp 


0.1 


0.999 ± 0.001 


0.98 ± 0.03 


1.00 ± 0.06 


1 


1.000 ± 0.001 


0.96 ± 0.02 


1.00 ± 0.04 


10 


0.999 ± 0.003 


1.00 ± 0.01 


0.95 ± 0.02 
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